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Abstract. Recently Hirota and Kimura presented a new discretization of the Euler 
top with several remarkable properties. In particular this discretization shares with the 
original continuous system the feature that it is an algebraically completely integrable bi- 
Hamiltonian system in three dimensions. The Hirota-Kimura discretization scheme turns 
out to be equivalent to an approach to numerical integration of quadratic vector fields 
that was introduced by Kahan, who applied it to the two-dimensional Lotka-Volterra 
system. 

The Euler top is naturally written in terms of the so (3) Lie-Poisson algebra. Here we 
consider algebraically integrable systems that are associated with pairs of Lie-Poisson al- 
gebras in three dimensions, as presented by Giimral and Nutku, and construct birational 
maps that discretize them according to the scheme of Kahan and Hirota-Kimura. We 
show that the maps thus obtained are also bi-Hamiltonian, with pairs of compatible Pois- 
son brackets that are one-parameter deformations of the original Lie-Poisson algebras, 
and hence they are completely integrable. For comparison, we also present analogous dis- 
cretizations for three bi-Hamiltonian systems that have a transcendental invariant, and 
finally we analyze all of the maps obtained from the viewpoint of Halburd's Diophantine 
integrability criterion. 

1. Introduction 

The problem of numerical integration, namely that of approximating the flow of a 
smooth vector field by an iterative scheme given in terms of a difference equation or a map, 
is one of the central problems of numerical analysis. If the underlying differential equation 
is Hamiltonian, or volume-preserving, or has some other important geometrical feature 
(such as being invariant under the action of a Lie group of symmetries), then as far as 
possible one would like to select a discretization scheme which preserves this feature, and 
this has led to the development of geometrical integration methods [4]. For the special case 
of completely integrable systems, ideally one would like to obtain discretizations which are 
themselves completely integrable. The area of integrable discretization has been developed 
quite extensively, especially from the Hamiltonian viewpoint, and a comprehensive review 
of the field can be found in the monograph [39]. In this paper we are concerned with a 
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novel approach to discretization, which was used by Hirota and Kimura to obtain new 
integrable discrete analogues of the Euler and Lagrange tops [15,22]. 

The discretization method studied in this paper seems to be introduced in the geometric 
integration literature by W. Kahan in the unpublished notes [20]. It is applicable to any 
system of ordinary differential equations for x : R — > IR n with a quadratic vector field 

x = Q(x) + Bx+ c, 

where each component of Q : W 1 — > W 1 is a quadratic form, while B e Mat nxn (lR) and 
c G R n . Kahan's discretization reads as 

^^ = Q(x,x) + ^B(x + x) + c, (1) 

where 

Q(x,x) = i[Q(x + x)-Q(x)-Q(x)], 

is the symmetric bilinear form corresponding to the quadratic form Q. Here and below 
we use the following notational convention which will allow us to omit a lot of indices: for 
a sequence x : Z — > R we write x for x^ and x for Xfc+i. Eq. (DO) is linear with respect to x 
and therefore defines a rational map x = /(x, e). Clearly, this map approximates the time- 
(2e)-shift along the solutions of the original differential system, so that x^ m x(2ke). (We 
have chosen a slightly unusual notation 2e for the time step, in order to avoid appearance 
of powers of 2 in numerous formulae; a more standard choice would lead to changing 
e i— > e/2 everywhere.) Since Eq. ([T]) remains invariant under the interchange x «-> x with 
the simultaneous sign inversion e i— > — e, one has the reversibility property / _1 (x, e) = 
/(x, — e). In particular, the map / is birational. 

Kahan applied this discretization scheme to the famous Lotka-Volterra system and 
showed that in this case it possesses a very remarkable non-spiralling property. Some 
further applications of this discretization have been explored in [21,36]. 

The next, even more intriguing, appearance of this discretization was in the two papers 
by R. Hirota and K. Kimura who (being apparently unaware of the work by Kahan) 
applied it to two famous integrable systems of classical mechanics, the Euler top and 
the Lagrange top [15,22]. Surprisingly, the Kahan-Hirota-Kimura discretization scheme 
produced integrable maps in both the Euler and the Lagrange cases of rigid body motion. 
Even more surprisingly, the mechanism which assures integrability in these two cases 
seems to be rather different from the majority of examples known in the area of integrable 
discretizations, and, more generally, integrable maps, cf. [39]. We shall use the term 
"Hirota-Kimura type discretization" for Kahan's discretization in the context of integrable 
systems. 

In the recent paper [32] the Hirota-Kimura integrability mechanism has been further 
investigated and its application to the integrable (six-dimensional) Clebsch system has 
been considered. The integrability of the Hirota-Kimura type discretization of the Clebsch 
system has been established, in the sense of: i) existence, for every initial point, of a four- 
dimensional pencil of quadrics containing the orbit of this point; ii) existence of four 
functionally independent integrals of motion. Note that for the purposes of paper [32], 
integrability of a dynamical system is synonymous with the existence of a sufficient number 
of functionally independent conserved quantities, or integrals of motion, that is, functions 
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constant along the orbits. Other aspects of the notion of integrability, such as Hamiltonian 
properties or explicit solutions, still require further investigation. However, it is known 
that algebraically completely integrable cases of geodesic flow on 5*0(4) are related to 
the intersection of four quadrics in P 6 [2]. The Hirota-Kimura method of discretization 
has been recently applied to the classical three-dimensional nonholonomic Suslov problem 
in [7]. 

The above examples of Hirota-Kimura type discretizations suggested the following 

Conjecture [32]. For any algebraically completely integrable system with a quadratic 
vector field, its Hirota-Kimura type discretization is algebraically completely integrable. 

Since algebraically completely integrable systems generically correspond to linear flows 
on abelian varieties [40] , this statement should be related to addition theorems for multi- 
dimensional theta-functions. 

The aim of this paper is both to study how this novel method of discretization applies 
to a set of algebraically integrable systems in three dimensions, and to see how these 
results compare with the analogous discretizations of some quadratic vector fields with 
transcendental invariants. The former set of systems considered are algebraically inte- 
grable in the sense that they have a sufficient number of algebraic integrals in involution; 
however, there are various other (more stringent) notions of algebraic complete integra- 
bility. Through this study we are able both to verify the above conjecture for a new set of 
examples, and to gain further understanding of how the integrability of the discretization 
depends on the algebraic nature (or otherwise) of the integrals of motion in the original 
continuous system. Kahan's discrete Lotka-Volterra system illustrates the subtlety of this 
dependence, as we now describe. 

Kahan used his approach to discretize the Lotka-Volterra system 

x = x(i-y), y = y(x-i), 

which preserves the Poisson bracket {x, y} = xy, or equivalently, the symplectic form 
l/(xy)dx A dy. This is an integrable system with one degree of freedom, x = {x,H}, 
y = {y, H} with Hamiltonian 

H = \ogxy — (x + y). 

Kahan's discretization for this system reads [20] 

x — x ^ _ _ y — y ~ ~ ~ /nS 
— - — = x + x - xy - xy, — - — = —y - y + xy + xy. (2) 

This discretization preserves the same symplectic structure as the original system of ordi- 
nary differential equations, for which it provides a numerically stable integration scheme 
which appears to retain the qualitative features of the continuous orbits (which are closed 
curves H = constant in the positive quadrant x > 0, y > 0) [37]. 

In fact, as noted in [31], Kahan's discrete Lotka-Volterra system is algebraically inte- 
grable for e = ±1. To be precise, when e = 1, it reduces to the second order recurrence 
x n +\Xn-i = x n (2 — x n ) for the x coordinate, which belongs to the class of antisymmetric 
QRT maps studied in [43]. This recurrence is linearizable (the iterates satisfy a linear 
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recurrence of sixth order [19]), and the map has the integral 

y z x z \x y J 

which (for fixed H) defines a quartic curve of genus zero; this is also an integral for e = — 1 
(which can be seen immediately from the reversibility property of Kalian's discretization 
scheme). However, for other non-zero values of e, Kahan's discrete Lotka-Volterra system 
should not be algebraically integrable; we present some numerical evidence for this in 
section 7 below. Indeed, since the integral H for the original system is transcendental, 
from continuity arguments one would expect that (at least for small enough e) any integral 
of the the discretization should be transcendental as well. Further numerical studies, as 
mentioned in [32], indicate that this discrete system may well be non-integrable, with 
characteristics of chaos only evident by zooming in deeply on regions of the phase plane. 

The outline of the paper is as follows. In the next section, we briefly review the Euler 
top together with the discrete Euler top found by Hirota and Kimura. In section 3 we 
describe six quadratic bi-Hamiltonian flows in three dimensions, which were presented 
in [11] (extending results in [3]), and are associated with pairs of real three-dimensional 
Lie algebras. Moreover, each of these systems, which we denote by Si for i — 1, . . . , 6, is 
algebraically integrable; the system £ 6 is equivalent to a special case of the Euler top. The 
fourth section is devoted to applying the Hirota-Kimura discretization scheme to these 
six systems, to obtain discrete systems (or maps) in three dimensions which we denote 
by dSi, and in section 5 we present the explicit solutions of these maps for i — 1, . . . , 5 
(the case i = 6 being already included in the work of Hirota and Kimura [15]). Section 6 
is concerned with applying the same discretization method to three other bi-Hamiltonian 
systems from [11] which have transcendental integrals. In section 7 we present the results 
of applying Halburd's Diophantine integrability test to each of the maps obtained, and 
prove that all but one of them are Diophantine integrable in the sense of [12]. The final 
section is devoted to some conclusions. 



2. Euler top and its Hirota-Kimura type discretization 

The 50 (3) Euler top is a well-known three-dimensional bi-Hamiltonian system belonging 
to the realm of classical mechanics [34]. The differential equations of motion of the Euler 
top read 

x = aiyz, y = a 2 zx, z = a^xy, (3) 

with «j being real parameters of the system. We recall that this system can be explicitly 
integrated in terms of elliptic functions, and admits two functionally independent integrals 
of motion. Indeed, a quadratic function H(x) = 71X 2 + ^y 2 + ls,z z is an integral for Eqs. 
([3]), if 7i«i +7202 + 72«2 = 0. In particular, the following three functions are integrals of 
motion: 

H% = a^y — 012Z , H 2 = a\z — a^x , H 3 = a 2 x — a\y . 
Clearly, only two of them are functionally independent because of a\H\ +a 2 H2 +03-^3 = 0. 
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The Hirota-Kimura discretization of the Euler top introduced in [15] reads as 

x — x = ea,\ (yz + yz) , 

y — y — ea 2 (zx + zx) , (4) 
J — z = ea 3 (xy + xy). 

Thus, the map / :xhx obtained by solving (Tj0) for x, is given by: 



x = /(x, e) =A x (x, e)x, A(x, e) = — ea 2 z 1 — ea 2 x . (5) 




Apart from the Lax representation which is still unknown, the discretization ([5]) exhibits 
all the usual features of an integrable map: an invariant volume form, a bi-Hamiltonian 
structure (that is, two compatible invariant Poisson structures), two functionally inde- 
pendent conserved quantities in involution, and solutions in terms of elliptic functions. 
For further details about the properties of this discretization we refer to [15] and [30]. 

3. Some bi-Hamiltonian flows related to real three-dimensional Lie 

algebras 

Hamiltonian systems in three dimensions provide the simplest non-trivial examples 
of degenerate Poisson structures, where the rank of the Poisson tensor is less than the 
dimension of the phase space. In three dimensions, a non-trivial Poisson tensor P has 
rank two at generic points of the phase space, which means that (at least locally) there 
exists a Casimir function K and another function <p such that 

{Xj,x k \ = e jH (p— (6) 

in local coordinates xi,x 2 ,Xs] cf. Theorem 5 in [9]. This can be expressed in invariant 
form, by using the standard volume three-form Q = dx\ A dx 2 A dx% in IR 3 to associate P 
with the one-form J = Pj Q = <p dK. An important thing to observe from the form of 
the Poisson bracket is that in three dimensions the Poisson tensor can be multiplied 
by an arbitrary function while preserving the Jacobi identity. 
Given an Hamiltonian system 

x={x,#} 

defined in terms of the bracket (J6]) with an Hamiltonian function H (functionally inde- 
pendent of K), it is clear that the equations of motion have two independent integrals, 
namely H and K. Moreover, by fixing the value of the Casimir function (which may not 
be defined everywhere), we can regard this locally as a system with one degree of free- 
dom which is integrable on each of the two-dimensional symplectic leaves K = constant. 
However, for complete integrability the global existence of H and K is required. 

Giimral and Nutku made a detailed study of the geometry of three-dimensional Pois- 
son structures, and considered the conditions for the existence of globally integrable bi- 
Hamiltonian structures [11]. For a given three-dimensional system to be bi-Hamiltonian 
it is necessary and sufficient that the Jacobian at an arbitrary point be a Poisson tensor 
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and that there exist two globally defined and (almost everywhere) functionally indepen- 
dent integrals of motion. Associated with two independent integrals K and H, there are 
two compatible Poisson tensors, such that K is the Casimir for one Poisson structure 
while H is the Casimir for the other. In other words, if the dimension is three then 
two compatible Poisson tensors are completely determined by the constants of motion, 
and according to a relevant Theorem by Magri [24], provided certain technical conditions 
are satisfied, bi-Hamiltonian systems are completely integrable in the sense of Liouville- 
Arnold. Furthermore, in this setting there is an invariant volume form Q (not necessarily 
canonical) which is preserved by the bi-Hamiltonian flow. An important example of such 
flows corresponds to Nambu mechanics [28], given by 

x = VH x VK, 

which in these coordinates gives a divergenceless vector field (divx = 0); this means that 
the canonical measure is preserved by the flow. In particular, the Euler top is an example 
of Nambu mechanics in three dimensions; for other examples of Nambu mechanics in 
optics and elsewhere, see [17]. 

In [11] the authors present a list of all non-trivial bi-Hamiltonian flows that are asso- 
ciated with pairs of real three-dimensional Lie algebras and their Casimir invariants (as 
described in [29]); this list extends results in [3]. To be more precise, they consider pairs 
of real Lie-Poisson algebras defined by pairs of linear Poisson structures P, Q, and write 
down vector fields x = V satisfying 

V = -PdK = --QdH, 

c 

where K is the Casimir for Q, while H is the Casimir for P (and minus signs are included 
in order to be consistent with the conventions of Giimral and Nutku). In [11] twelve such 
systems are presented, extending a list in [3], and each flow preserves a corresponding 
measure given in coordinates (xi, X2, £3) = (x,y,z) by 

Q = c dx A dy A dz, 

related to the standard volume form by the conformal factor (or multiplier) cJ3 

To begin with, we shall be concerned with only six out of the twelve systems in Giimral 
and Nutku's list, namely the ones which have non-transcendental integrals of motion. 



In fact, on page 5704 of [11] the authors state that the given systems are all "with multiplier unity", 
and denoting the multiplier by M they say "these equations have M = 1 [...] they are Nambu mechanics 
representatives" , but as should be clear from Table 1 this is not the case: two of the systems given there 
have a non-constant multiplier. For systems with non-constant multiplier, it is remarked in [11] that 
they may be only locally (but not globally) equivalent to Nambu mechanics, by a suitable change of 
coordinates. 
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They read 



Si : 


x = 


-x 2 , 


y = -xy, 


z 


= 2y 2 + xz ; 


(7) 


£2: 


x = 


-x 2 , 


y = xy, 


z 


= — 2y 2 + xz ; 


(8) 


£3: 


X = 




y = -yz, 


z 


2 1 2 
= x +y ; 


(9) 


£4: 


X = 


—xz, 


y = yz, 


z 


2 2 

= x -y ; 


(10) 


£5: 


X = 


xy, 


y = -x 2 , 


z 


= y(2x - z) ; 


(11) 


£ 6 : 


X = 


y(2z - x), 


y=x - z , 


z 


= y(z - 2x) . 


(12) 



These six systems are examples of algebraically completely integrable systems, in the sense 
that in each case the integrals are algebraic (in fact, rational) functions of the coordinates 
x, y, z for which the Poisson structures P, Q are linear. (We consider two other examples 
where there is one transcendental invariant in section 6.) 

The corresponding linear Poisson structures and integrals of motion are given in Table 
[TJ For instance, the flow Si, given in Eq. (J7J), admits the bi-Hamiltonian structure given 
by the compatible pair (P^\ c^Q^ 1 '), where 

P« : = {x, y} = 0, P$ = {y, z} = y, P« = {z, x} = -x, 

: = {x, y} = x, Q$> = {y, z} = z, Q$ = {z, x} = 2y, 

with conformal factor c\ = 1/x 2 . The quantities Hi = y/x and Ki = zx + y 2 , preserved 
by the flow, are respectively the Casimir functions of P^ and This is equivalent to 
say that the following Lenard-Magri chain [24] is satisfied: 

P^dHi = 0, P^dKi = -Q^dHi = -(-x 2 , -xy, 2y 2 + xz) T Q^dK x = 0, 

Cl 

where dHi and dKi denote the differentials of the functions Hi and Ki respectively. The 
same scheme holds for the flows Si with 2 < i < 6. 

In Table [1] there are actually just five independent Lie-Poisson structures, namely 

p(l) = p(3) ; p(2) = p(4) ; p(5) 5 p(6) = g(l) = q(2) = g(5) ; g(3) = g(4) = g(6) > CQrre _ 

sponding respectively to the Casimir functions Hi = H 3 , H 2 = H4, H 5 , H 6 = K\ = K 2 = 
K 5 , K 3 = = Kq. The last column in Table [T] gives the associated real three-dimensional 
Lie algebras; see [29] for more details. Observe that the flows S\ and Sq each correspond to 
a particular case of the equations of motion ([3]) of the so (3) Euler top. More precisely, for 
S4 one has to make the change of variables (x, y,z) t— > (x — y, x + y, z) and fix the param- 
eters so that (aii, «2, 03) = (— 1, — 1, 1), while for £ 6 one takes (x, y,z)\-^(x — z,x + z, y) 
and the parameters are (oti, a 2 , a 3 ) = (1, —3, 1). 

4. HlROTA-KlMURA TYPE DISCRETIZATION OF THE FLOWS £j 

The goal of this section is to show that Hirota-Kimura type discretizations of the bi- 
Hamiltonian flows Si, 1 < i < 6, provide completely integrable discrete-time systems. 
The following result holds. 

Theorem 1. The Hirota-Kimura type discretizations of the bilinear flows Si, 1 < i < 6, 
given in Eqs. ^ffUfy. read 

d£i- x = i4 j - 1 (x;e)x = A(x;-e)x, (13) 
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Table 1 . Lie-Poisson structures, invariants, conformal factors and related 
real three-dimensional Lie algebras 



i 


pi 
*12 




pi 


Wl2 


^23 


W31 


TJ 


Ki 


Ci 


9 


1 





y 


—x 


X 


z 


2y 


V 


zx -\- y 2 




1 


A 3 , 3 , *I(2,R) 
















X 






X 2 


2 





-y 


—x 


X 


z 


2y 


xy 


zx + y 




1 


e(l,l),5t(2,R) 


3 





y 


—x 


z 


X 


y 


y 


\(x 2 +y 2 + 


z 2 ) 


1 


^3,3, so(3) 
















X 




X 2 


4 





-y 


—x 


z 


X 


y 


xy 


\(x 2 +y 2 + 


z 2 ) 


1 


e(l,l),so(3) 


5 





X 


y 


X 


z 


2y 


\{x 2 +y 2 ) 


zx + y 2 




-1 


e(2), S [(2,R) 


6 


X 


z 


2y 


z 


X 


y 


zx + y 2 


\{x 2 +y 2 + 


z 2 ) 


-1 


sl(2,R), so(3) 



where the matrices A(x;e) are given in Tabled The quantities Hi(e),Ki(e), given in 
TablelE, are integrals of motion for the maps ( T73l) . Moreover the maps (T73|) preserve the 
volume form 



a 



H 4 K, 



■dx Ady A dz, 1 < i < 6. 



(14) 



Note that for small e the birational maps (Tl3l) approximate the time shift along the tra- 
jectories of the corresponding continuous equations of motion (|71TT2l). The same invariant 
volume form fTl4l) . which is independent of e, is preserved by both the continuous and the 
discrete systems. 

Proof: We shall prove Theorem 1 for just one case, namely i = 5. The remaining 
cases can be proved by similar straightforward computations. 

The Hirota-Kimura discretization of the flow £5, given by Eq. (fTTl) . reads explicitly as 



x — x = e(xy + xy), 
y — y = —2exx, 

z — z = e{2xy + 2xy — yz — yz) , 



(15) 



that is 
with 

A 5 (x; e) 

The fact that the quantities 

H 5 (e) - 



A 



5 x (x; e)x = A 5 (x; -e)x, 



1 — ey —ex 
2ex 1 

— 2ey —e(2x — z) 1 + ey 



1 x 2 + y 2 



21 + e 2 x 



2ry,2 ' 



KM 



zx + y 



1 + e 2 x 



2^2 ' 
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v4j(x; e) 

'l + 2ex 
ey 1 + ex 
—ez — 4ey 1 — ex , 

1 + 2ex 
— ey 1 — ex 
— ez 4ey 1 — ex . 

'1 +ez 


-2ex 



1 + 


-2ex 



1 — ey —ex 
2ex 1 

—ley —e(2x — z) 1 + ey / 

1 + ey e(x — 2z) —2ey 

-2ex 1 2ez 

2ey e(2x — z) 1 — ey. 






ex\ 


1 + ez 




-2ey 


V 





ex ' 


— ez 




2ey 





#1 



#5 



1 — e 2 x 2 



1 + e 2 (x 2 + y 2 ^ 



1 + e 2 x 2 



H 6 
1 — 3e 2 xz 



e 2 x 2 



K, 



1 — e 2 x 2 



l + e 2 (x 2 + y 2 ^ 



1 + e 2 (x 2 + |/ 2 ) 



K R 



1 + e 2 x 2 



^6 

1 — 3e 2 xz 



are integrals of motion of the map (j!5p is proved by the following computation. Equation 
i?s(e) = H 5 (e) means that 

{x-x){x + x) + (y-y)(y + y) = -e 2 (xy -xy)(xy + xy), 

that is, using Eq. ( floT) . 



-(xy + xy)(x + x) — xx{y + y) — —-{xy — xy){x — x), 

which is an algebraic identity. A similar computation shows that equation K 5 (e) = K 5 (e) 
is identically satisfied. 

We now prove that the map (fl5l) preserves the volume form 



(xz + y 2 )(x 2 + y 2 ) 



dx A dy A ob, 
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which is equivalent to saying that 

83c _ (xz + y 2 ) (x 2 + y 2 ) 
<9x (xz + y 2 )(x 2 + y 2 ) 

First of all we note that differentiating Eq. ( Tl5l) with respect to x, y, z one obtains the 
columns of the matrix equation 

e) — = A>(x; -e). 

Computing determinants lead to 

<9x det A 5 (x; -e) (1 - ey) (1 + ey + 2e 2 x 2 ) 



det 



<9x detA 5 (x;e) (1 + ey) (1 - ey + 2e 2 x 2 ) ' 

Now, by using the map (Tl5l) . a straightforward computation shows that the relation 

(xz + ^(x 2 + y 2 ) _ (1 - ey) (1 + ey + 2e 2 x 2 ) 
(xz + y 2 )(x 2 + y 2 ) (1 + ey) (1 — ey + 2e 2 x 2 ) 

holds identically. 

□ 

In the construction of an invariant Poisson structure for the maps (TT~3T) we shall make use 
of results from [5] (Proposition 15 and Corollary 16 there), which we restate here. Suppose 
that / : M — > M is a smooth mapping of an n-dimensional manifold M, with an invariant 
volume form Q (that is, f*Q = Q). Define to to be the dual n- vector field to Q such that 
ujjQ = 1, where as usual the symbol j denotes the contraction between multivector fields 
and forms. It follows that if Ii, . . . , I n -2 are integrals of / with dI\A- ■ -A/ n _2 7^ 0, then the 
bivector field a = u_i dl\ ■ ■ -j rf/ n _ 2 is an invariant Poisson structure for /. If J 1? . . . , J„_ 2 
is another set of independent integrals and r = uja dJ x ■ ■ • j <iJ„_2 is the corresponding 
Poisson structure, then a and r are compatible, i.e. for any constants a, b, the bivector 
field acr + br is again a Poisson structure. 

In particular for n — 3, if a three-form Q, given by Eq. (114n in our case, is invariant 
under a map / defined by (EE]), we can define the dual trivector field 

.Odd 
ox ay oz 
so that for any integral I of / the bivector field 

JT .(did d dl d d dl d d 

a = uj^dl = (p{x, y, z) I A — + — — A — + — — A — 

\ oz ox Oy ox Oy oz oy oz ox 

is an invariant Poisson structure for /, as well as any linear combination of such bivector 
fields. Explicitly, the Poisson brackets of coordinate functions are given by 

{x, y} = 0(x, y, z)^-, {y, z] = <f>{x, y, z)^-, {z, x} = <f>{x, y, z)^j-. 

Note that the inverse volume density 4>(x, y, z) can be multiplied by an arbitrary integral 
of / without violating the Poisson property. 
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For the maps ( fl3l) . the invariant Poisson structures P^\ c t can be computed 
according to the following formulae: 



1 HK 



d HA e) KAe 



' Ejkl 



dHjje) 
dxi 



and 



— Qk(o = - 



HK 



dKAe) 



d HA e) KAe 



£jkl 



dxi 



(16) 



(17) 



with 1 < i < 6, the summation convention is assumed for the index /, and above we have 
used (xi,x 2 ,x 3 ) to denote (x,y,z). (The reader should note that these indices 1,2,3 for 
the coordinates in K 3 should not be confused with the index n used to denote iterates of 
maps in subsequent sections.) This corresponds to taking = HiKi/ci above, and then 
rescaling by the inverse of the product of the integrals, l/Hi(e)Ki(e), in each case. Thus 
the following statement holds. 

Theorem 2. The maps (T73|) admit the compatible pair of invariant Poisson structures 
(P«(e),c- 1 Q«(e)) ; where P«( e ) and Q«(e) are given respectively in Tables and \^ 
The conformal factors q are the same as in the continuous case, given in Table d 



Table 3. First deformed Poisson structure 



1W) 











2 2\ 

e x ) 



2„,2\ 



x(l + 3e z y 



y(i 

-y(l +e 2 x 2 ) 
y[l + e 2 (x 2 + y 2 )} 
-y[l-e 2 {x 2 ~y 2 )] 

x (1 — e 2 y 2 ) 

z(l + 3eV) 



2 9> 



—x (1 — e x 
—x (1 — e 2 x 2 ) 
-x[l + e 2 (x 2 + y 2 )} 
-x [1 + e 2 (x 2 - y 2 )} 
y(l + e 2 x 2 ) 
2y(l-3e 2 xz) 



Note that Eqs. (I16H17P provide one-parameter deformations of the Lie-Poisson tensors 
given in Table [U This is equivalent to saying that Tables [3] and H] provide 
deformations of the real three-dimensional Lie algebras y4 3)3 ,s[(2,R),so(3), e(l, 1), e(2). 
Finally we note that the integrable discrete-time system dS^ is just a particular case 
of the Hirota-Kimura discretization of the so (3) Euler top [15], whose bi-Hamiltonian 
structure has been presented recently in [30]. 
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Table 4. Second deformed Poisson structure 



i 




Q23(e) 


Qai(e) 


1 


X 


z + e 2 x(^x + 2y 2 ) 
1 — e 2 x 2 


2w 


2 


/ -i 9 9 \ 

x (1 — e z x^) 


9 / r\ 9 \ 

z + e z x(;za; + 2|/ z ) 


2y (1 - e 2 x 2 ) 


3 


z 


x (1 — e^z) 


ft 9 \ 

1/(1 ~e 2 z) 


l +e 2( x 2 +2/ 2) 


l + e 2 (x 2 + y 2 ) 


4 


z[l + e 2 (x 2 + y 2 )} 


x(l-e 2 z 2 ) 


y(l-eV) 


5 


x (1 + e 2 x 2 ) 


z — e 2 x(xz + 2y 2 ) 


2y(l + eV) 


6 


z + §e 2 x(x 2 + y 2 — 2 2 ) 


x + |e 2 2;(2; 2 + y 2 — x 2 ) 


y (1 — 3e 2 xz) 



5. Explicit solutions to the integrable systems dS i7 1 < i < 5 

As shown in [15,22], and recently in [30], the integrable discrete-time systems obtained 
through the Hirota-Kimura type discretization seem to admit a straightforward construc- 
tion of their explicit solutions, at least for the case of three-dimensional maps. Here we 
provide the explicit solutions for the discrete-time integrable systems dSi with 1 < i < 5. 
The cases i = 4, 6 are each special cases of the 50 (3) Euler top, whose solutions, both con- 
tinuous and discrete, are investigated in [15,30], so here we present the solution only for 
% = 4, since i = 6 is similar. For comparison, in Table [5] we give the explicit solutions for 
the continuous-time flows £i with 1 < i < 5. The parameters a, (3, 7, 9, A, /i, k appearing 
in the table can easily be expressed in terms of the initial conditions and/or the integrals 
of motion. 

We now construct the explicit solutions to the discrete-time systems dSi with 1 < i < 5, 
thus providing the discrete counterpart of Table El Let us recall that we consider each of 
x, y, z as functions on eZ. To simplify the notation we set x = x n , y = y n , z = z n , so that 
x = x n+ i, y = y n +i, z = z n+ \. For the sake of brevity, henceforth the discrete integrals of 
motion Hi(e) and Ki(e) in Table [5] will be denoted respectively by Hi and Ki, 1 < % < 6. 

The following statement holds. 



Theorem 3. The explicit solutions to the integrable maps dSi, 1 < i < 5, given by Eq. 
(E^, read: 



THREE-DIMENSIONAL DISCRETE SYSTEMS OF HIROTA-KIMURA TYPE 

Table 5. Solutions to continuous systems S iy 1 < i < 5 



13 



i 


x{t) 


y(t) 


z(t) 


Integrals 


I 


1 


a 






( H 1= f3 


t + a 


t + a 




\ K l=1 


2 


1 

t + a 


[3 [t + a) 


(t + a) 






( H 2 = (3 

1 7V 

A2 = 7 


3 


\cos9 


\sa\9 


Atanh[A(t + a)} 




( H 3 = tan 6> 


cosh[A(i + a)] 


cosh[A(i + a)] 




1 ^3 = A 2 /2 


4 


fdn [X(t + a)] 
+fcn [X(t + a)} 


fdn \X(t + a)] 
-■fen [\{t + a)\ 


fcAsn [X(t + a)] 


r # 4 = a 2 (i - k 2 )/4 

\ K 4 = A 2 (l + k 2 )/4 


5 


Asech [X(t + a)\ 


-Atanh[A(i + a)] 


/iCosh[A(t + a)] 
+Asech [X(t + a)] 


/ H 5 = A 2 /2 
\ X 5 = A( M + A) 



• i = I: 



1 



Vn = 



2e(n + r)' 

P 

2e(n + r)' 



z n = 27e (n + r) - 



P'e- 1 + 7 6 
2(n + r) ' 



with Hi — (5, Ki — 7; 
i = 2: 



n 2e(n + r)' 



(n + r) 



2(n + r) 



4(n + r) 
1 

2(n + r) 



[ 7 -/3V(4(n + r) 2 -l)], 



wit/i #2 = (3, K 2 = 7; 
i = 3: 



cos # sinh <5 



x 



ecosh(25n + k) ' 
sin 6? sinh <5 
^™ ecosh(25n + n) ' 
z n = e _1 tanh 5 tanh(25n + k) , 



(18) 
(19) 
(20) 



(21) 
(22) 
(23) 



(24) 

(25) 
(26) 
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with 



H-i = tan. 6, 



tanh 5 
2e 2 ' 



4: 



sn 5 
~2~T 

sn 6 



dn(2n5 + k) k cn(2n5 + k) 



2e 



cn 5 dn S 

dn(2n5 + k) kcn(2n5 + n 



cn 5 



dn 5 



sn 5 sn(2n5 + «) 



where sn, cn, dn are i7ie Jacobian elliptic functions with modulus k [42], and 
(l-k 2 )sn 2 5 - 1/1 cn 2 o"dn 2 o" 



#4 



5: 



2e 2 [2- (l + P)sn 2 o"] 

sinh 5 
ecosh(2o"ra + k) ' 



2& 



2 cn 2 5 + dn 2 o" 



with 



—e 1 tanho~tanh(2o~n + k), 

\i cosh(2o~n + k) + (e _1 + /isinh<5) sinhtfsech {25n + k), 



— — tanh 2 5, K 5 = e 1 // sinh 5 + e 2 tanh 2 S. 
2e 2 



(27) 

(28) 
(29) 

(30) 

(31) 

(32) 
(33) 



Proof: Let us illustrate the procedure to find the solutions (I18H33P for just one of the 
five discrete systems dSi, 1 < % < 5. We shall consider % = 3. The remaining cases can 
be verified by elementary direct computations, apart from dS^, which we reserve for the 
Appendix. 

The system d£ 3 reads: 



Un+1 ~ Un = ~^(Un+l z n + Vn z n+l)i 
Z n+1 ~ z n — 2e(x n +iX n + UnVn+l)- 



(34) 
(35) 
(36) 



It has two integrals of motion, 



Vn+l 



l + e 2 (x 2 n + y 2 n 



X n+1 + Un+1 + Z n+1 

l + e 2 (x 2 n+1 +y 2 +1 ) 



in involution with respect to the pair (P^(e), l Q^(e)), as given in Tables [3] and HI 

The solution to the continuous-time flow £ 3 , as in Table El suggests the following ansatz 
for the solution of the map 



ucosO 
cosh T„ 



z/sin# 
cosh T„ 



A tanhT n , 



(37) 
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with constant parameters A, u, 9. By substituting the ansatz (1371) into the formulae for 
the integrals, we see that H 3 = tan#, while 

^ _ A 2 + (y 2 - A 2 )sech 2 T ra 
3 2(1 + e 2 i/ 2 sech 2 T n ) 

is constant (for all T n ) if and only if v 2 = A 2 / (1 — e 2 A 2 ). Upon setting A = e _1 tanh5, in 
terms of another parameter 5 (with 8/e — 0(1) in the continuum limit e — > 0) this gives 
v 2 = e~ 2 sinh 2 5 and K3 = sinh 2 5/(2e 2 ). Substituting the ansatz into the third part of 
the map, namely ( l36i) . and using the addition formulae for hyperbolic functions, one can 
see that this equation implies that 

sinh(T n+ i - T n ) = ^-^r / = sinh 25, 
tanho 

hence T n+1 — T n = 25. This implies that T n = 25n + k for some constant k, and then it 
is straightforward to verify that Eqs. (1341) and (1351) are also satisfied identically. 

□ 



6. Discretization of three-dimensional bi-Hamiltonian flows with one 

transcendental invariant 

There have been several studies of integrable Hamiltonian systems which have tran- 
scendental invariants [8,13]. Among the six bi-Hamiltonian flows with transcendental 
invariants listed in [11] we select the following ones: 



£7: 


X = 


-x 2 , 


V = 




z = 2£y 2 + xz ; 


(38) 




X = 


-x\ 


V = 


-x(x + y), 


z = 2y[x + y) + xz ; 


(39) 


£ 9 : 


X = 


—xz, 


V = 




z = x 2 + £y 2 . 


(40) 



In [29] the real parameter £ is restricted to the range |£| G (0, 1), but here we need not 
impose this requirement. Observe that the equations of motion ( 138]) reduce to the flow 
£\ if £ = 1 and the flow £ 2 if £ = — 1. Also, the equations (140]) reduce to £ 3 if £ = 1, and 
to £ 4 if £ = -1. 

The Lenard-Magri chains for the flows (I38H40|) are given by 

P^dHi = 0, P^dKi = -Q®dHi Q (i) dKi = 0, 

Ci 

for % = 7, 8, 9 respectively, with = = and = (related to sl(2, R) 
and to so (3) respectively, see Tabled]), 

p(?) : P $ = fa y} = o, p(J) = fa z} = fa P g) = fa x} = - Xj 

P®: P®=fay} = 0, P£ ) = {y,z} = x + y, = {z, x} = -x, 
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and P (9) = P (7) , with 



H 7 


= yx ? , 


K 7 


= H G 


Hs 


= xe~ y/x , 


K 8 


= H 6 


H 9 


= yx't, 


K 9 


= K 6 



zx + y 2 , cj = x 

-y/x 



.2 



e 



zx + y , c 8 - 

x 

^-(x' + y' + z 2 ), c 9 = x-^. 

Thus the transcendental invariants in each case are given by H7, H 8 and H 9 respectively. 
(Strictly speaking, H 7 = H 9 is only transcendental when £ (jL Q, otherwise it is algebraic.) 
Moreover, note that the Lie algebra related to is actually a one-parameter family of 
Lie algebras, parametrized by £; see [29] for more details. It can also be regarded as a 
four-dimensional Lie algebra, by taking y = log y as a new coordinate and regarding £ as 
a central element. 

We now construct the Hirota-Kimura type discretizations of the flows £7, £$ and £9; 
these are denoted using the notation introduced in section [5l 

6.1. Explicit solutions to d£j. The explicit solution to the equations of motion ( 1381) is 
given by: 

x(t) = -l—, y(t)=P(t+a)-e, z(t) = (t + a)[ 1 -(3 2 (t + a)-^], (41) 
t + a 

with H7 = (5 and K7 = 7. Following the approach described in section [51 the discrete-time 
version of the flow £7 reads: 

x n j r \ x n l€x n x n j r \ , (42) 

Vn+l —y n = -^(Un+lXn + y n X n +l), (43) 

z n +i ~ z n = e(x n+1 z n + x n z n+1 ) + Ae£y n y n+ i. (44) 

The decoupled equation for x n can be rewritten as a total difference, 

1 1 

= 2e, 



Xn+l X n 

from which it follows by summation that 

Xn = , V' T = O ' 45 

2e(n + r) 2ex 

this is the discrete version of the first equation in (1411) . to which it tends in the continuum 
limit 

e — > 0, 2m — > t, 2er — > a. 

By substituting x n given by Eq. ( )45l) into Eq. ( 1431) we get a difference equation for the 
variable y n , whose solution reads 

r(3 r(n + l + r-£/2)r(r + £/2) 
yn 2e(n + r)r(n + r + £/2)r(r + l-£/2)' 1 } 

where T(z) is the complete gamma function. We can now solve Eq. ( H6l) for the constant j3 
(up to scale) to write it as a function of x n and y n , which gives an explicit transcendental 
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integral: 

y n T^/2 + (2ex n )- 1 ] 



H 7 



x„r(l-e/2 + (2ex n )- 1 )' 
Now inserting x n and y n , given respectively by Eqs. (I45H46I) into Eq. (jHJ) we find a 
difference equation for z n . Its solution is 

r 7 [4(n + r) 2 -l] 2r 2 /^[4(n + r) 2 - l]r 2 (r + g/2) 
n 2e(4r 2 - l)(n + r) e r 2 (r + 1 - £/2) (n + r) n ' 1 J 

where 

™ _ y [2( J + l + r)-e]r 2 (j + l + r-e/2) 

" PO' + T ) + 3 ] (20' + r ) + t 4 0' + ^ - !] F2 0' + r + e/2) ' 

In principle, Eq. (1471) can implicitly be solved for 7 (after first replacing j + r by (2ea;.,) _1 
everywhere to remove explicit dependence on the parameter r), to give another transcen- 
dental invariant K 7 , in that case the bi-Hamiltonian structure can be reconstructed by 
the same formulae as above in cases 1-6; this means that the system d£? is completely 
integrable. Using the formula for H-j above we can reconstruct one invariant Poisson 
bracket for this map explicitly, as 

{x,y} = 0, 

- i/2 + (2ex)- 1 ) - v&(£/2 + (2ex)- 1 ] 



{y,z} = 2/(1 - eV) 
{z, x} = x(l — e 2 x 2 ) 



2ex 



1 



where \I/ is the digamma function. This bracket has Hj as a Casimir, and for £ = ±1 (up 
to scaling) it reduces to the brackets P"(e) and P^ 2 \e) respectively. 

It is straightforward to verify that the explicit form of the solution for x n ,y n , z n given 
by Eqs. (I45H47I) can be used to recover the previous formulae for the discrete systems dS\ 
and d?2 given in Eqs. ( fl8ll20l) and ( [2TH23]) by setting £ = ±1 in the respective cases. 



6.2. Explicit solutions to d£$. The explicit solution to the equations of motion (159]) is 
given by 

X(t) = -—, y(t) = ^—± L, z (t )= 7 t+a - l ^—± ^, 

t + a t + a t + a 

with H 8 = e _/3 and K 8 = 7. The discrete-time version of the flow £ 8 reads: 

x n+i ~~ Xn — — 2ex n x n +i, (48) 
y n+ i -y n = -e(y n+1 x n + y n x n +i) - 2ex n x n+1 , (49) 
z n +i ~ z n = e(x n+1 z n + x n z n+ i) + 2e(y n+1 x n + y n x n+1 ) + Aey n y n+l . (50) 

The first equation for x n is identical to that in the previous case, and has the solution 
x n = (n + r) _1 /(2e) as before. By substituting x n into Eq. (|49|) we get a difference 
equation for the variable y n , whose solution reads 

Vn = 777 — : — v , (51) 
2e(n + r) 
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where 

U n = (n + t + 1/2) - * (r + 1/2) , 

with \I/(z) denoting the digamma function as before. This leads to the transcendental 
invariant 

H a = ^ + Hf {1/2 + (2ex n )~ 1 ). 

Upon inserting x n as in (j4"51) and y n given by Eq. (15T1) into Eq. (150]) we find a difference 
equation for z n , whose solution is given by 

t [4n(n + 2r) (/3 2 r + /? + 7) + 4r 2 7 - 7] 2 [4(n + r) 2 - 1] . 



2e (4r 2 - 1) (n + r) 2e(n + r) 



-K , (52) 



where 



n-l 



v 1 + 2/3r + ^[2(j + r)(l + 2gr) - 1 + 2gr] - f/ 2 [2(j + r) + 1] 
" ^ [2(j + r)+3][2(j + r) + l][4(j + r)2-l] 

Similarly to the situation for <i£ 7 , the system d£$ has another transcendental integral 
Kg which is given implicitly by solving Eq. ( 1521) for 7. This implies that <i<?8 is also 
bi-Hamiltonian and hence completely integrable. 

6.3. The system d£g. For all values of the parameter £, the equations of motion ( HOl) 
can be reduced to a quadrature, namely 

r x (t) d s 
t + const = ± / — ; . 

7 sV2H -s 2 - K 2 s 2 ^ 

Given x(t) determined by this quadrature, y and z are then given by 



y(t) = Kx(ty, z(t) = ±y/2H - x(t) 2 - K 2 x(t) 2 t . 

The constants H and K are respectively the values of H$ and Kg along an orbit. For 
certain values of £ the quadrature can be performed explicitly; for instance, when £ = 1 
it becomes an elementary integral, and the problem reduces to the solution of £3, while 
for when £ = — 1 it becomes an elliptic integral, corresponding to the solution of £4, as 
given in Table 5. The case £ = 1/2 is also an elementary one, while £ = —1/2 and £ = 2 
also give elliptic integrals (of the first and third kind, respectively). More generally, for 
all rational values of £ this quadrature is an hyperelliptic integral. 

However, it is straightforward to check that the cases £ = ±1, which were solved 
already, are the only ones for which the system has the Painleve property (i.e. all solutions 
are meromorphic functions of t in these cases only). In general the solutions have movable 
algebraic branch points in the complex t plane when £ G Q, and movable logarithmic 
branch points when £ ^ Q. 

The qualitative nature of the solutions is fairly insensitive to the parameter £. In fact, 
for £ > the trajectories interpolate between the two fixed points (x, y, z) = (0, 0, ±y/2H), 
at the north/south poles of the sphere x 2 + y 2 + z 2 = 2H, while for £ < there are closed 
periodic orbits. These two types of behaviour are exemplified by each of the explicitly 
solvable cases £ = ±1. 
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The Kahan-Hirota-Kimura discretization of this flow is given by 

Vn+1 ~ Vn 
2-ra+l Zji 

We have not attempted to solve this discrete system in the case £ ^ ±1. In fact, numerical 
results for the latter case (as described in the next section) provide evidence for the non- 
integrability of the system for generic values of £. 

7. DlOPHANTINE INTEGRABILITY TEST 

Over the past fifteen years or so there has been a gradual development of methods for 
testing integrability of maps or difference equations, using such concepts as singularity 
confinement [10], algebraic entropy [14], Nevanlinna theory [1] and orbit counting over 
finite fields [35]. In certain limited cases it has been proved that these tests provide 
necessary conditions for integrability of a map, in a suitable sense, most usually in the 
setting of algebraic integrability (see [23], for instance), but in general it is an open 
problem to determine when these tests are effective. 

Most recently Halburd proposed an extremely simple criterion for integrability which 
applies to rational maps defined over Q (or more generally over a number field), which 
he named the Diophantine integrability test [12]. For a map whose n-th iterate has 
components x n G Q, written as a fraction x n = p n /c[n in lowest terms, the height of x n 
is defined to be H(x n ) = max(|p n |, |g n |); this is the archimidean height of x n , and the 
logarithmic height is h(x n ) = \ogH(x n ). For a map in dimension N, with N components, 
the height H n of the n-th point on an orbit is defined to be the maximum of the heights 
of all the components, with h n = logH n being the logarithmic height. Halburd defined a 
map to be Diophantine integrable if the logarithmic height h n of the iterates of all orbits 
has at most polynomial growth in n. If we define the Diophantine entropy along an orbit 
O to be 

E(0) := lim -log/i n , 

n^oo n 

then a Diophantine integrable map is one for which E(0) = for all orbits. 

Diophantine entropy is somewhat similar to algebraic entropy [14], which measures the 
height growth of rational functions generated by rational maps. In the latter setting the 
height of each iterate is just the maximum of the degrees of the polynomials in the nu- 
merator and denominator, considered as a rational function in the initial data. However, 
a huge disadvantage of using algebraic entropy is that one must usually try to guess a 
recursive relation to generate the degrees of these polynomials. The great advantage of 
Halburd's test is that it is extremely quick and straightforward to implemement numer- 
ically with a computer, and if the map is Diophantine integrable then a plot of log/t n 
against logn should look asymptotically like a straight line (see Figured]), otherwise it 
will have an exponential shape (see Figure CJ). The main drawback of using the test is that 
at present it has the status of a distinct definition of integrability, and it is not clear how 
it is related to other such definitions, like complete integrability in the Liouville- Arnold 
sense. 



€yX n Z n ^_i -\- X n +{Z n ) , 
= -^(Un+lZn + y n Zn+l), 

= 2ex n+l x n + 2e£y n+1 y n . 
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FIGURE 1. Plot of log h(x n ) versus logn for the first 200 iterates of the 
integrable map d£ 3 for e = 1/2 with initial conditions x = 3/7, y = 
11/13,20 = 23/47. 

Despite these drawbacks, it is worth remarking that, at least for maps in two or three 
dimensions, Diophantine integrability is a necessary condition for algebraic integrability. 
For example, a two-dimensional map which is algebraically integrable has a conserved 
quantity whose level sets are algebraic curves. Assuming that each of these curves is 
irreducible, and that not all orbits of the map are periodic, it was observed by Veselov [41] 
that they must all have genus zero or one; this follows from a theorem of Hurwitz which 
says that curves of genus two or more have automorphism groups of finite order [27]. (This 
argument also extends to the case when the level curves are reducible.) If the curve is 
rational (genus zero), then the map can be linearized, in which case the logarithmic heights 
grow linearly, h n ~ Cn for some constant C, while a curve of genus one is birationally 
equivalent to an elliptic curve, for which the heights grow as h n ~ Cn 2 . (See chapter 17 
in [6] for an introduction to archimidean heights on elliptic curves, or chapter VIII in [38] 
for a more general discussion of heights.) Similar considerations apply to algebraically 
integrable maps in three dimensions, where the algebraic curves are the level sets of two 
independent integrals, or to systems with N — 1 algebraic integrals in N dimensions (as 
considered in [23] from the viewpoint of singularity confinement). However, in general 
these level sets can have two or more irreducible components; see [16] for several examples 
with two components in three dimensions. 

Here we prove that all of the discrete systems constructed here, except for dSg, pass 
the Diophantine integrability test, before presenting numerical results which show more 
detailed behaviour of the growth of heights for some of these systems. For the theoretical 
and numerical analysis here it is convenient to set e = 1/2; since the right hand sides of 
the difference equations are homogeneous (of degree two), this can always be achieved by 
scaling x n , y n , z n by the same factor. 

Theorem 4. The discrete systems dSi for i — 1, ... ,8 are all Diophantine integrable. 

Proof: Without loss of generality we set e = 1/2, as mentioned above, and consider 
each of the maps with rational initial data xo,yo,%o (and parameter £ e Q for the case 
of d£i). This implies that all of the iterates (x n ,y n ,z n ) of these birational maps are 
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also rational numbers for all n (except on a set of initial data where these maps become 
singular) . 

For the maps d£\ and d£<i it is clear from the explicit solutions, as given in Eqs. (TTS1 - 
|2"U1) and Eqs. (I2HI23P respectively, that in each case the iterates are given in terms of 
parameters a, [3, 7 G Q, and these rational iterates have numerators and denominators 
which grow linearly in n. Hence the logarithmic height satisfies h n = logn + 0(1) (sub- 
polynomial growth) for these two maps. 

The maps d£% and d£$ are naturally considered together, because their explicit solutions 
given in Theorem [3] are in terms of hyperbolic functions (or equivalently, exponential 
functions of n) in each case, which means that the intersections of the level sets of their 
two integrals are curves of genus zero. This implies that the heights of iterates should 
grow like h n ~ Cn. To prove this directly for g?£ 3 , note that one can eliminate y n from 
Eq. ( I36i) by setting y n = H 3 x n , and then further eliminate z n between that equation and 
Eq. ( 1341) to get an expression of the form z n = F(x n , x n +i, H 3 ) with F being a rational 
function. This leads to a single recurrence of second order for w n = l/x n , namely 

4w 3 n+1 + (1 + H%)(2w n+1 + w n ) 

W n+ 2 = — • 

4w n w n+1 - (1 + Hi) 
The latter recurrence has the conserved quantity 

1 _ 2{wl + wl +l ) + l+^Hl 

Aw n w n+1 - (1 + Hi) 

and furthermore admits the linearization 

w n+2 - 2Lw n+1 + w n = 0, (53) 

which linearizes the system d£%\ in terms of the original integrals and solution parameters 
we find L = (2 + K^)/{2 — K%) = cosh 25. From the second order linear recurrence 
(|53|) it follows directly that the height H(w n ) grows exponentially with n, and hence 
h{w n ) = h(x n ) ~ Cn (cf. Figure [T]) for some C > 0. Since y n = H 3 x n , and z n can be 
written as a rational function of x n and x n+ i, it follows that h(y n ) and h(z n ) also have 
linear growth in n. Analogous arguments apply to d£$. 

Similarly, it is natural to consider the maps d£± and d£ & together, because the intersec- 
tions of the level sets of their two integrals are curves of genus one; the details for d£± are 
given in the Appendix. For d£± each of the coordinates x n ,y n ,z n of a point on an orbit 
can be written in terms of Jacobi functions, which are related by a Mobius transformation 
to the Weierstrass p function. For instance, the solution for z n in (1291) is linear in the 
Jacobi sine, which is an elliptic function of order two with two simple poles in each period 
parallelogram; this implies that a relation of the form z n = {aX n + b)/ (cX n + d) holds, for 
some constants a, b, c, d, where X n is the nth term in a sequence of X coordinates of points 
Po + nP ^ E, for an elliptic curve E given in Weierstrass form as Y 2 = X 3 + AX + B (for 
some A, B). It is known that, as long as P is not a torsion point (which would correspond 
to a periodic orbit), the height grows like h(X n ) ~ Cn 2 as n — > 00, where the constant 
C > only depends on the height of the point P [38]. Since z n is related to X n by a 
rational map of degree one, it follows that h(z n ) has the same quadratic growth in n, and 
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FIGURE 2. Plot of log h(x n ) versus logn for the first 125 iterates of the 
integrable map d£^ for e = 1/2 with initial conditions x = 7/3,y = 
11/13, z = 23/47. 



similarly for h(x n ) and h(y n ). The same arguments apply to d£$, this being a special case 
of the Hirota-Kimura discrete Euler top, whose solutions are most naturally written in 
terms of Jacobi functions. 




Figure 3. Plot of log h(y n ) (bottom set of points) and log h(z n ) (top set 
of points) versus logn for the first 2000 iterates of the integrable map d£ 7 
for e = 1/2 with initial conditions x = 3/7,y = 11/13, z = 23/47 and 
parameter!; = 19/17. The bottom points have been fitted against log n+ 1.64 
(a straight line on this scale), and the top points against logn + 2.64; the 
curve logn + log logn is also shown. 



Finally, for the systems d£ 7 and d£$ we make use of direct estimates of the growth of 
heights, based on the original maps. For both these systems, note that from the explicit 
solution we have h(x n ) = h{n + r) = logn + 0(1). It is convenient to define Y n = y n /x n 
and Z n = z n /x n in each case, and then note that h(y n ) = h(Y n ) + O(logn), and similarly 
for h(z n ). From the second part of the map d£ 7 we have 
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which implies 

h(Y n+1 )-h(Y n )<logn + 0(l) =^ h{Y n ) < n\ogn + 0(n), 

where the second implication follows by summing over n. Thus h(Y n ) has weaker than 
quadratic growth in n. Similarly for Z n we have 

(n + r- l/2)Z n+1 = {n + r + 3/2)Z n + 2£Y n Y n+1 , 

which implies that 

h(Z n+1 ) < h{Z n ) + h(Y n Y n+1 ) + O(logn) < h(Z n ) + nlogn + 0(n), 

and hence h(Z n ) < logri + 0(n 2 ), which is weaker than cubic in n. For d£$, analogous 
estimates show that h(Y n ) < nlogn + 0(n) and h(Z n ) < 2n 2 log n + 0(n 2 ), so this system 
is Diophantine integrable as well. 

□ 

Having proved that the systems are all Diophantine integrable, we can compare the the- 
oretical results with some numerical experiments. For the system d£$ we see that the 
log- log plot gives what we expect: genus zero means linear growth of logarithmic height, 
so \ogh(x n ) = logra + 0(1); this is evident from the plot of points in Figured], which lie 
asymptotically on a straight line of slope 1. Similarly for the genus one case, we expect 
\ogh(x n ) = 21ogn + 0(1), and Figure [T] shows points which asymptote to a line with 
slope 2. In this case the offset, corresponding to the correction at 0(1), is function of the 
height of a point on an associated elliptic curve, and both the point and the curve vary 
with the initial data of the map. 




2.86- „ 

2.84- 

2.82- 



Figure 4. Plot of log(h(x n )/n) versus n for d£^ with the same data as 
Figure El 

The theoretical results on the growth of heights for the algebraically integrable sys- 
tems dSi for 1 < i < 6, as detailed in the above proof, are confirmed by the numerical 
calculations, and for those cases we have an exact expression for the leading order as- 
ymptotic behaviour. Moreover, one can also look at how the asymptote is approached. 
Taking the system dE^ for example, h(x n )/n approaches a constant as n —>■ oo, and from 
the numerical plot in Figure H] one can see that this limit is reached in a very uniform 
manner, in keeping with a correction of 0(l/n) to this constant. Similarly, in the case of 
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dS^, corresponding to motion on an elliptic curve, we see from Figure [5] that once again 
the convergence of h(x n )/n 2 to a constant appears to be almost monotone. 




Figure 5. Plot of log(h(x n )/n 2 ) versus n for 0IS4 with the same data as 
Figure 



The non-algebraically integrable cases, dSj and d£g, have some extremely interesting 
features compared with the others. First of all, the method of proof used in Theorem 
4 above has not necessarily provided the leading order asymptotics of the logarithmic 
heights, but has merely given upper bounds on the growth of the form Cn J log n with 
j — 1 for h(y n ) and j = 2 for h(z n ) in each case. Let us focus on the case of dS^. 
Upon looking more closely at Eq. (j54p . it would appear that the upper bound for h(Y n ) 
might be sharp, so that h(Y n ) ~ nlogn (and h(y n ) would have the same leading order 
asymptotics). However, studies of particular sequences of rational iterates show that 
cancellations occur between the numerator and denominator of Y n and the prefactor 
(n + a + 1 — £/2)/(ra + a + £/2), which means that the height of Y n+ \ is therefore smaller 
than the crudest estimate for the upper bound. This weaker growth has a knock-on 
effect, meaning that the growth of heights of z n also seems to be much weaker than 
expected. Indeed, Figure E] suggests that the correct asymptotics should be linear growth 
in n for the logarithmic heights of both y n and z n , i.e. h(y n ) ~ C\n and h(z n ) ~ C^n for 
positive constants C^Cg. For the particular sequence of heights plotted in that figure, 
a numerical fit shows that logCi ~ 1.64 and logC*2 ~ 2.64 rs logCi + 1. We have also 
plotted log n + log log n for comparison, to show how the upper bound for h(y n ) fails to be 
sharp. Another surprising feature of this system is that for different choices of initial data 
we find (to within numerical accuracy) the same values of C\ and C2', this is in contrast to 
the algebraically integrable setting described above, where the coefficient in front of the 
leading order term is dependent on the initial data. Thus we might conjecture that for 
this map C\, C2 are independent of initial data, and also that C2 = e C% holds identically, 
in which case there should be some deeper arithmetical explanation for this asymptotic 
behaviour. Similarly to the case of d£f, numerical results for the system dSg also show 
linear growth of logarithmic heights. 

Supposing that the numerical observation of linear growth of h n for these discrete sys- 
tems with transcendental invariants is indeed the correct asymptotic behaviour, it is then 
interesting to look at how h n /n approaches a constant value. The results we find are in 
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Figure 6. Plot of 400 iterates of log(h(x n )/n) versus n for dEj with the 
same initial data and parameters as in Figure 



stark contrast to the algebraic setting: rather than the almost monotone convergence seen 
in the previous examples, for d£ 7 we find that h n /n shows rapid fluctuations which persist 
for increasing values of n. These fluctuations in the asymptotics are somewhat reminiscent 
of the "random" -looking error terms that appear in some famous arithmetical functions, 
such as the difference between the prime-counting function n{n) and the logarithmic in- 
tegral [25]. It would be interesting to know whether these fluctuations might provide a 
means of characterizing the difference between discrete systems which are algebraically 
integrable and those with transcendental invariants. 




Figure 7. Plot of \ogh(x n ) versus logn for the first 17 iterates of Ka- 
han's discretization of the Lotka-Volterra system with initial conditions 
xo = 1/2, yo = 5/3 and parameter e = 9/14. 

For comparison with the Diophantine integrable examples above, in Figure [7J above we 
have plotted the growth of h n for a particular case of the discrete Lotka-Volterra system 
due to Kahan, which is the degree two birational map given in Eq. (EJ). This figure shows 
that the logarithmic height seems to grow exponentially, indicating non-integrability of 
this system. Indeed, the heights of iterates grow so fast that even on a fairly new computer 
it took 1 hour to calculate the heights of 17 rational iterates with Maple; the value of 
h{xn) is of the order of 1Q 325009 in this case. Upon examining the data used in Figure [7] 
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more carefully, it is apparent that h(x n+ i) ~ 2 h(x n ) to very good accuracy, so we expect 
h n ~ C 2 n . This would mean that the logarithmic height essentially doubles with each 
step, giving a Diophantine entropy of log 2 for generic (aperiodic) orbits. This appears to 
be the same as the algebraic entropy of the map when e 2 7^ 1, which was calculated by A. 
Ramani [33]. In any case, this is the entropy value that one would expect for a generic 
(non-integrable) birational map of degree two. 

Finally we should mention the results of numerical calculations of the growth of heights 
for the map dSg for various cases with £ 7^ ±1 (that is, excluding the two special cases 
where the map is already known to be algebraically integrable) . For generic rational values 
of the parameter £ we find that the map dSg is not Diophantine integrable, but rather the 
Diophantine entropy is log 3 for generic orbits, this being the typical value to be expected 
for a non-integrable birational map of degree three. (See Figure [H] for an illustration of 
the exponential growth of logarithmc heights when £ = —1/2.) These numerical results 
suggest that while that the original continuous system £$ is algebraically integrable (in the 
sense of having two independent algebraic integrals), the corresponding discrete system 
is not. We shall return to this point in our conclusions. 




Figure 8. Plot oflogh(x n ) versus logn for the first 11 iterates of the 
map dSg for e = 1/2 with initial conditions xq = 2/5, yo = 7/3, = 11/13 
and parameter £ = —1/2. 



8. Concluding remarks 

In this paper we studied three-dimensional birational maps which provide integrable 
time-discretizations of quadratic bi-Hamiltonian flows associated with pairs of real three- 
dimensional Lie algebras, as presented in [3, 11]. We have shown that for the six cases of 
continuous flows which are algebraically integrable, the Hirota-Kimura type discretization 
provides maps admitting two independent rational integrals of motion, in involution with 
respect to a pair of compatible Poisson tensors. We have also provided explicit solutions 
of the resulting discrete systems dSi for i = 1, . . . , 6, which are given in terms of either 
rational, hyperbolic or elliptic functions in each case. These results confirm the conjecture 
that the property of algebraic integrability is preserved by this discretization scheme. 

We have also applied the same procedure to two cases of integrable continuous flows in 
three dimensions having one rational and one transcendental integral of motion, for which 
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the resulting maps, d£ 7 and d£ 8 , admit explicit solutions in terms of rational functions 
and either gamma or digamma functions. In each of these cases, we have found an explicit 
formula for one transcendental integral of the map, but the second integral is only defined 
implicitly by the solution. Nevertheless, this is sufficient to assert that the latter two cases 
are also completely integrable in the Liouville- Arnold sense. Therefore the Kahan-Hirota- 
Kimura discretization scheme preserves integrability even in these transcendental cases. 
However, for another example of a continuous integrable system with one transcendental 
integral, it appears likely that the corresponding discretization d£$ is not integrable for 
generic values of the parameter in the map. 

In an attempt to gain a better understanding of the difference between the algebraic 
and transcendental cases, we have also analyzed all of these discrete systems from a 
different viewpoint, within the arithmetical setting of Diophantine integrability. So far, 
the Diophantine integrability test has been applied to various algebraically integrable 
systems and discrete Painleve equations [12], as well as to certain birational maps that 
are not algebraically integrable and fail the test [18,19]. Our theoretical results show 
that d£i and d£$ provide examples of discrete integrable systems with transcendental 
invariants that are also Diophantine integrable, in the sense defined by Halburd. Moreover, 
more detailed numerical results suggest that these discrete integrable systems might be 
distinguished from algebraically integrable maps by the manner in which the logarithmic 
heights converge to their leading order asymptotics. This asymptotic behaviour deserves 
to be studied more carefully in the future. 

A further interesting point is the comparison with Kahan's discretization of the Lotka- 
Volterra system, which is an integrable flow in the plane with a transcendental integral. 
The discrete system provides a non-standard symplectic integrator of this flow, and seems 
to preserve the qualitative features of the continuous counterpart. However, the numerical 
results indicate that the discrete system is not Diophantine integrable, which adds further 
evidence to the conjecture that (for generic values of e) it should not be Liouville integrable 
either. 

Similar considerations apply to the discrete system d£ 9 : numerically it appears that 
for generic values of the parameter £, it fails the Diophantine integrability test. Since 
the continuous system £g has algebraic integrals for ( 6 Q, this means that in general 
algebraic integrability (in the weakest sense of the term) is not preserved by the Kahan- 
Hirota-Kimura discretization. Actually one could already observe this for the system £7, 
since it has the integral H-j which is algebraic for all £ G Q, but the integral H 7 for d£j is 
transcendental unless £ G Z. This suggests that one should impose a much stronger notion 
of algebraic complete integrability (a.c.i.) if this is to be preserved by the discretization 
scheme. For instance, one can require that the generic level sets of the integrals are 
smooth abelian varieties, possibly extended by (C*) m for some m. An excellent discussion 
of various different definitions of a.c.i. can be found in chapter V of [40]. 

As for the case of the discrete three-dimensional Euler top (see [15,30]), there is one 
other standard attribute of integrable systems that remains to be found for the maps d£i 
for i — 1, ... ,8, namely their Lax representation. This is an open problem which deserves 
further investigation. 
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There are three other bi-Hamiltonian flows in the list of Giimral and Nutku, all of which 
have one transcendental invariant. Preliminary results suggest that their Kahan-Hirota- 
Kimura type discretizations are qualitatively similar to the system d£ 9 , and we expect 
that these maps are not Liouville integrable. We reserve the study of these systems for 
future work. 

Finally we remark that another non-standard symplectic integrator for the Lotka- 
Volterra model, with similar numerical properties, has been given by Mickens in [26] . It 
would be very interesting to see if the approach to discretization proposed by Mickens 
shares some of the remarkable properties of Kahan's. 
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Appendix: solution of the system d£ 4 in elliptic functions 

Here we derive the formulae ( 12711291) corresponding to the solution of dS^, as in Theorem 
EJ The first observation to make is that the curve V in affine three-dimensional space 
defined by the equations 

V: xy = H [l + e 2 (x 2 + y 2 )] , x 2 + y 2 + z 2 = 2K [l + e 2 (x 2 + y 2 )] , 

corresponding to the intersection of the level sets H± = H, K4 = K, has genus one (at 
least for generic values of H and K). To see this, note that V is a double cover of the 
curve 

C: xy = H[l + e 2 (x 2 + y 2 )] 
in two dimensions, via the covering map 

tt: V -> C 

(x,y,z) i-> (x,y) 

which is ramified over the four points (x, y) G C obtained from the simultaneous solutions 
of xy = H(l + e 2 (x 2 + y 2 )), (1 - 2Ke 2 )(x 2 + y 2 ) = 2K (when z = 0). Since the curve C is a 
conic (genus zero), it follows from the Riemann-Hurwitz formula [27] that V is (the affine 
part of) a curve of genus one. The first order recurrence relations for x n ,y n ,z n , namely 

•En+l ^(•^n+l^n ~\~ ^n^n+l)) 

Un+l ~ Un = t{yn+l z n + Un z n+l) ? 
z n+l ~ z n = 2t( x n+l x n ~~ Z/n2/ra+l)) 

correspond to a birational map from this curve to itself, inducing an automorphism of an 
isomorphic elliptic curve (i.e. a plane curve defined by a Weierstrass cubic), and it follows 
that x n = X(u + nv ) for a suitable elliptic function X, and similarly for y n , z n . One can 
see some points on a real connected component of such a curve in Figure [91 



(55) 
(56) 
(57) 
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FIGURE 9. Plot of the first 1000 points on the orbit of the integrable map 
0IS4 with initial conditions Xq = 7/3, yo = 11/13, Zq = 23/47. 

From the equations for V it is easy to see that the functions corresponding to x n ,y n , z n 
each have simple poles at the same places, and they are elliptic functions of order two. 
These facts suggest that it may be most convenient to write the formulae in terms of 
Jacobian (rather than Weierstrassian) elliptic functions. Indeed, if we set 

s = x + y, d = x — y, 

then the equations for V become 

(1 - 2He 2 )s 2 - (1 + 2He 2 )d 2 = AH, (1 - 2Ke 2 )(s 2 + d 2 ) + Az 2 = AK, 

which are reminiscent of (linear combinations of) the quadratic relations 

sn 2 (w) + cn 2 (w) = 1, k 2 sn 2 (u) + dn 2 (u) = 1, 

for Jacobi functions. 

In order to obtain the formulae ( 127H291) . it is instructive to take a detour through Jacobi 
theta functions, by deriving bilinear equations from Eqs. ( 15511571) . Upon setting 

A n + B n ^ A n — B n C n 

2eD n 2eD n eD n 

the system ( !55ll is equivalent to the following three bilinear equations: 

A n +iD n — A n D n+ \ = —(B n+ iC n + B n C n+ i), 
B n +iD n — B n D n+ i = — (A n+ iC n + A n C n+ i), (58) 
Cn+i-Dn — C n D n+ \ = A n+ iB n + A n B n+ \. 

One should hesitate to call (1581) the Hirota bilinear form of Eqs. ( I55ll57j) . because there 
are four unknowns (tau-functions) A n , B n ,C n , D n but only three equations, so the system 
is underdetermined. Despite this apparent problem, we can solve this bilinear system in 
terms of Jacobi theta functions j = 1, ... ,4, by comparing these equations with the 
identities in exercise number 3 on page 488 of [42]; the first of these is the relation 

± v)$ 2 (u =F v)^ 4 = tf^M^M^Mv) ± Mu)Mu)M v )Mv) (59) 
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the last is 

tf 3 (w ± v)$ A {u =F v)tf 3 #4 = Mu)Mu)M v )Mv) T &x(u)Mu)M v )M v ) ( 60 ) 

and there are four other relations of this kind, for different permutations of the four indices. 
Here $j without argument denotes a theta constant (i.e. i?j = $j(0), which depends on the 
modulus k). By taking the sum of the two equations given in with opposite choices of 
± signs, and similarly taking the difference of the two equations specified by ( |60l) . one sees 
that both ■di(u + v) , &2(u — v)+~&i(u — v)'&2{u+v) and ^(u+v^^u — v) — #3 (« — ?;) $4 («+?;) 
are proportional to i?i(w)#2(w), modulo -^-dependent factors. Thus if v is regarded as a 
fixed constant, and the shift u — > -u + 2v is identified with n — > n+1, then the first equation 
in ( 1581) is satisfied if (suppressing all arguments and the index n) the identifications 

C~^l, fi~^2, ^4-^3, £>~#4 

are made, up to suitable f-dependent scaling denoted by the ~ symbol. Moreover, these 
identifications are consistent with the second and third equations in (|58p . which are conse- 
quences of the aforementioned other four bilinear relations between Jacobi theta functions. 
Given that the Jacobian elliptic functions are defined in terms of theta functions by 



AM^iV ^Mu/^IV ^MW^iV 

it follows that the solution of the difference equations (1531) has the form 

x n = A dn(« + 2n5) + \i cn(« + 2nd), 
y n = A dn(/t + 2nS) — fi cn(/t + 2nS) , 
z n = z/sn(/€ + 2nS), 

for constants 5, k and suitable prefactors A, /i, z/ which are given in terms of 6 and the 
theta constants. The expressions ( 12711291) can also be verified directly from the addition 
formula for sn, namely 

. , sn(u)cn(t>)dn(t>) + sn(v)cn(u)dn(n) 

sn(w + t;) = , (61) 

1 — k z sn z {u)sn z {v ) 

as well as the analogous formulae for cn and dn. Using ( 16TT) to calculate sn(u+v )— sn(w— v), 
and then setting w — > k+ (2n + 1)5, t> — > 5, gives an expression for the left hand side of the 
third equation in (1551) . and performing analogous computations for the right hand side 
and for the other two difference equations allows the prefactors A,/i, v to be determined 
directly in terms of Jacobi functions with argument 5, in agreement with (I27H29I) . 

Finally, note that the solution depends on the required number of arbitrary constants, 
namely the three parameters 5, k, k. The parameter 5 and the modulus k are determined 
by the values of the integrals if 4 and K±, by solving the relations fl30l) as a system for k 
and sn<5 and then performing the elliptic integral 5 = f^ 5 d^/rj, while k is found from 

K= /-a* „ 2 =(i-a(i-m 

Jo V 
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